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Abstract. A set of equations is derived from the Boltzmann kinetic equation de- 
scribing charge transport in semiconductors. The unknowns of these equations de- 
pend on the space-time coordinates and the electron energy. The non-parabolic 
and parabolic band approximation are treated in detail. In these cases, the set 
of equations is equivalent to those obtained in the spherical-harmonics expansion. 
Stationary and homogeneous solutions are explicitly treated. In order to solve nu- 
merically the equations, a cut in the energy range is introduced. The modified model 
maintains the physical characteristics of the original equations. The solution of an 
asymptotic equation is found and compared to the numerical solutions. 
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1. Introduction 

The motion of electrons in a crystal is governed by complex physical laws so that 
accurate models are required to achieve correct results. The drift-diffusion model 
has been widely used in the past, and it has provided a good description of relevant 
physical mechanisms. In modern devices, whose sizes are in the submicron range, 
non-equilibrium effects play an important role and are not adequately modeled by 
the drift-diffusion approach. Since in many processes a more accurate description 
than the hydrodynamical setting is required, Boltzmann equation or Monte Carlo 
simulations are employed. A fully kinetic treatment of carrier dynamics guarantees 
accurate results but requires very expensive numerical procedures in order to solve 
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realistic problems. 

To reduce the complexity of the use of the full Boltzmann equation, many authors 
(see ref. and references therein) have introduced simpler models, assuming par- 
ticular forms for the probability density function. The aim is to yield results, which 
are more accurate than those obtained by hydrodynamical models but less expensive 
than direct particle simulations or numerical treatment of the Boltzmann equation. 
In this framework, a well-known model is derived using a spherical-harmonics ex- 
pansion 0, IQ, p from the Boltzmann equation. 

The paper is organized as follows. In Sec. 2, we introduce the Boltzmann Equa- 
tion. The collisional operator takes into account the interactions between electrons 
and phonons, i.e. vibrations of the lattice. This is assumed to be in thermal equilib- 
rium. A set of kinetic equations is derived in Sec. 3. These are obtained for a general 
expression of the microscopic kinetic energy, so that, for example, they can be used 
also for non-parabolic bands. The equations are of partial differential and difference 
type. Equations are specialized in Sec. 4 in the parabolic band approximation, and 
is Sec. 5, where the Kane band model is assumed. In both cases we show that the 
equations are equivalent to those obtained in the spherical harmonics framework. 
In Sec. 6 we write the equations in the stationary and homogeneous case. The rest 
of the paper is devoted to the study of these equations. In Sec. 7 an asymptotic 
equation, valid for large value of the electric field, is found and analyzed; the explicit 
solution is found. In Sec. 8 we discuss some difficulties related to the infinite energy 
range and we propose a suitable solution by introducing a cut in the scattering rate 
of the collisional operator. The numerical scheme to solve the equations is presented 
in Sec. 9, and the results are shown in the last section. 

Throughout the paper, boldface and lightface symbols denote vectors and scalar 
quantities respectively. 
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2. Basic equations 

We consider an electron gas, which moves in a lattice, subjects to an external electric 
field E. This can be applied or related, by Poisson equation, to the density of the 
gas. For the electron gas in a semiconductor device the Boltzmann equation |Q 
writes 

df 1 e 

(1) ^ + ^Vk^ • V./ - -E ■ Vk/ = Qif), 

where the unknown function /(t,x, k) represents the probability of finding an elec- 
tron at the position x G M^, with the wave- vector k G M^, at time t. The parameters 
h and e are the Planck constant divided by 2tt and the positive electric charge, re- 
spectively. The symbol Vk stands for the gradient with respect to the variables k 
and Vx that with respect to the space coordinates x. The particle energy e is an 
assigned nonnegative continuous function, defined on M^, of the wave- vector k. Here 
we assume that the low-density approximation holds, so that Q is linear in /. If 
K(k, k) is the sum of the scattering kernels, which describe the nature of the inelas- 
tic collisions (for example, electron-polar or non-polar optical phonon scattering), 
and ii'o(k, k) is the kernel for the elastic collisions (for example, electron-impurity 
scattering), then the coUisional operator writes 0, 0, 

(2) Qif) = (n, + 1) / Kik, k)5(£(k) - e(k) - huj)fit, x, k) dk 

+ n„ / K(k, k) 6(e(k) - e(k) + nuj)f(t, x, k) dk 

+ f Ko(k, k) S{e{k) - £(k))/(t, X, k) dk - u{k)f{t, X, k) , 

where 

(3) z/(k) = n„ / K{k, k) 6{6{k) - e{k) - Huj) dk 

+ (n„ + 1) / K{k, k) S{e{k) - e{k) + hu) dk+ f Ko{k, k) S{e{k) - s{k)) dk. 

JR^ JR^ 
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Often, in the following we omit to write explicitly that / depends also on t and x, 
whereas we do not operate with these variables. 

The constant quantity rig represents the thermal-equilibrium number of phonons and 
is given by 

1 

" (^) - 1 

where ks is the Boltzmann constant and Tl is the constant lattice temperature. 

The symbol 6{e(k)—e(k)±huj) means the composition of the real- valued function 
£(k) — e(k) ± huj and the Dirac distribution 6. Since the function e may have 
many different expressions, it is not possible in general, by using a unique explicit 
technique, to transform each of the integral operators (§) in equivalent operators 
without Dirac distributions. Therefore, we shall maintain the original form of (|^) in 
the following. 

As usual, due to the isotropy of phase-space, the kernels ^(k, k) and -K'o(k, k) 
are symmetric with respect to the wave-vectors k and k. Moreover, since the kernels 
are scalar, they depend on k and k only through the scalar quantities ^(k), £:(k) and 
k • k. This implies in particular that the collision frequency u depends on k only 
through the variable e. 

3. Energy- kinetic equations 

In the framework of the Boltzmann equation for a perfect gas, a classical procedure, 
due to Grad 0], to reduce the dimension of the space of the independent coordinates, 
is the moment method. It consists in expanding the ratio between the probability 
density / and a local maxwellian function. Inserting this expansion in the Boltzmann 
equation, an infinite sequence of differential equations are obtained. In general 
no finite set of equations is uncoupled from the rest of the system. To obtain a 
determined system up to certain order m (usually 13) only m equation are considered 
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and the coefficients of tlie expansion of higher order are assumed neghgible. 

In order to obtain a new set of equations, where the unknown functions depend 
on the space-time coordinates and the microscopic kinetic energy, we follow the 
scheme analyzed in |^. This scheme partially recalls the moment method. 

If J'(k) is a polynomial and u a real variable, we multiply the Boltzmann equation 
(|l]) by !P(k)(5(£:(k) — u) and then formally integrate with respect to the variable k 
over the whole space M'^. It is possible to show |^ that the following equation can 
be derived: 

(4) -/^^/(t,x,k)T(k)5(e(k)-«)dk 

+ Vx / ^ v(k)/(t, X, k)T(k) 5{e{k) - u) dk 
-cE • — / fit, X, k)y(k)v(k)5(£(k) - u) dk 

+|e ■ fit, X, k) [Vky(k)] 5(e(k) -u)dk = j^^ g(/)T(k) 5(e(k) -u)dk. 
Here the molecular velocity is defined by means of the formula 

(5) v(k):=lVke(k). 

For any choice of the function T(k) a new equation is obtained. As it happens for 
the moment method, each equation contains, except for particular cases 0, at 
least two unknown functions. For example, if we choose T(k) = 1, then the left 
hand side of eq. contains in general the following unknowns 

(6) iV(t,x,M):= ( fit,^,k)5ieik)-u)dk, 

(7) V(t,x,M):= / /(t,x,k)v(k)5(£(k)-M)rfk. 

The scalar quantity Nit, x, u) is the probability density function to find an elec- 
tron with energy u at time t and position x. The function is non-negative. 



because the probability density / is always non-negative. We note that the pres- 
ence of the delta-function implies that and V vanish for all {t, x, u) such that 
u ^ Uq = min |£:(k) : k G M^j. If the function is known, then we can obtain 
all the hydrodynamical scalars. For example, to find the hydro dynamical density 
p(t,x), it is sufficient to integrate with respect to all values of energy u; i.e. 

r+oo 

p(t,x):= / N(t,:s.,u) du. 

The physical meaning of the other quantity V(t,x, m) can be understood with a 
similar argument. 

In the rest of this article we assume the function e(k) to be even [i.e. e(k) = 
£(— k)) and the kernels K and Kq to depend only on e(k) and ^(k). The first 
hypothesis is verified for the most common models used in the numerical simulations. 
The second are introduced in order to simplify the treatment of the equations. Let 

%{e{k),e{k)) := K{k,k) and Xo{e{k),e{k)) := Ko{k,k) . 

Then eqs. ®-(|D become 

(8) Q(/) = (n, + l)X{e{k),e{k) + huj)N{t, x, ^(k) + hu) 

+ nq%{eik) - nu,e{k))N{t,^,e{}i) - hu) 

+ %o{e{k),e{k))N{t, X, ^(k)) - p{e{k))f{t, x, k) 

with 

(9) z/(£(k)) = (rig + l)X{e{k) - hu, e{k))a{e{k) - hu) 

+ ng3<;(£(k), £(k) + hij)(r{e(k) + huo) + %Q{e{k),e{k))a{e{k)) , 

where 

(10) a{u) := f S{e{k)-u)dk 
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is the density of states. 

We are interested in deducing two equations in the only variables and V. 
By choosing J'(k) = 1 and J'(k) = v(k) eqs. (|^), (^) and (|TDp give the following 
equations 



dN dY 



dV ^ 



-v{u)V, = 1,2,3). 



where 

(13) G{N) := (n, + l)%{u, u + huj)a{u)N{t, x, m + hu) 
+ nq%{u, u — huj)(T{u)N{t, — huj) 

— [ng%{u, u + ?iuj)a{u + Tiu) + (n^ + l)%{u, u — liuj)a{u — huj)] N(t, x, u) . 

(14) u^u) = nq%{u,u + hu)a{u + hu) 

+ (rig + l)X{u, u — Tiuj)a{u — Tiuj) + %q{u, u)a{u) , 

(15) n,, := / f(t,x,k).;,(k)t;,(k)5(£(k)-«)rfk, 

dv^ik) 



(16) A,,:= / /(t,x,k)^^5(£(k)-n)rfk. 



''J 

Also the tensor Ajj is symmetric, due to eq. (^. We note that the effects of the 
elastic collisions, through the kernel 3Co, disappear in the operator G{N). 
The set of eqs. (p!T|)-(|T^) contains and V but also two tensors Aij and Iljj. There- 
fore, except in the case of spatially homogeneous solutions with null electric field, we 
need to assume some relations (usually called constitutive equations in the thermo- 
dynamic framework), which link Aij and Uij with A^ and V. These relations should 
also depend on the form of the microscopic kinetic energy e. 

We feel that at this stage it would be worthwhile to make a simple observation. 



We note that the hnearity of the Boltzmann equation is maintained in eqs. (^TJ)- 
(|T2|). To maintain this feature, we suggest to assume constitutive relations of the 
kind 

(17) Ilij{t,yi,u)=Pij{u)N(t,:s.,u) + Pi{u)Vj(t,:>i,u) + Pj{u)Vi(t,yi,u) 

(18) ^ij{t, X, u) = dij{u)N{t, X, u) + di{u)Vj(t, x, u) + dj{u)Vi(t, x, u) , 

where the symmetric tensors pij and dij and the vectors pi and di must be deter- 
mined. In the next sections we propose a simple choice in the cases of the parabolic 
and non-parabolic band approximation. 

4. The parabolic case 



The simplest and widely used expression for the microscopic kinetic energy JIT 



IS 



(19) ^("'^^^^ 

where m* is the value of the effective electron mass in the parabolic mass approxi- 
mation. In this case the tensor Aij reduces to 

(20) A,, = ^iV(t,x,M)5,, (^,j = 1,2,3), 

where 6ij is the Kroneker symbol; so that no further assumption on eq. (plSf ) is 
required. The tensor Uij is determined by assuming that ([TTD holds with Pi{u) = 
and moreover that the constitutive relation is exact if f depends on k only through 
the variable e. For /(t,x, k) = /(t, x, £:(k)), we have 

iV(t,X,M) =/(t,x,M)cr(M) 

nij(t,x,M) = /(t,x,M) / fi(k)t;j(k)5(e(k)-M)rfk (i,j = 1,2,3), 



so that 



a{u)pij{u) = I fi(k)fj(k)(5(e(k) — m) (ik. 



Now, it is easy to verify that 

where 6 is the Heaviside step function. Then, for every u > 0, 
f 2 u 

/ , ^i(k)t^j(k)5(£:(k) -u)dk = -—S^Ja{u) = 1,2,3). 

Hence, eq. (|T2p becomes 

r B 1 



d 

2uV^N - 2tE—(uN) + 3tEN 

ou 



dt 3m* 

In the following the functions and V will be the unknowns, which must satisfy 
eqs. ([Tl|) and (|2l|) . It is evident that the knowledge of these quantities cannot provide 
the probability density / uniquely. Also, it is very hard to prove that, if and V, 
at the initial time, derive from a probability density /, there exists a non-negative 
probability density, which gives A^ and V for all times. 

Actually the situation is similar to that of the moments in the hydrodynamical 
equations derived by the Grad method. Here usually, not even for the initial data, 
the existence of a probability density, which gives all the moments occurring in the 
equations, is considered. 

Therefore, we limit to assuming only the simple conditions that A^ is nonnegative 
and both the variables A^ and V vanish as m < 0. 

We note that these unknowns are simply related to the first two terms of the 
spherical harmonic expansion of the distribution function 

(22) fit, X, k) = /o(t, X, £(k)) + k ■ fi(t, X, £(k)) + ... . 

It is a simple matter to show that 

Ar(t,x,£(k)) 
/o(t,x,.(k)) = ^^-^, 

Mt,x,.(k))= ^^^^^^^^^^^^^ V(t,x,.(k)). 
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The corresponding equations for /o and fi (see ref. [Q) coincide with our eqs. (pH]) 
and (0). This is a consequence of the choice on the approximation of the tensor 
Uij. It is evident that different choices give different sets of equations. In this paper 
we do not investigate this possibility. 

Similar equations were proposed by Hansch in ref. [13| in the framework of quan- 



tum kinetic transport. Also in that case a problem of equating the number of equa- 
tions and unknowns arose. It was solved by simple physical arguments, but the 
equations are different from eqs. ([TTD -(pll). 



5. The non-parabolic band approximation 

For large values of electron energy eq. ([T9|) is not adequate. The following suitable 
expression was given by Kane 0, [pJ]) : 

(23) e{l + a,e) = ^k\ 

where ak is a constant. Now, it is simple to verify that 

/ / — r\ 3 



So eq. (p!T| ) can be written explicitly. The second equation requires only to determine 
Uij and Aij according to eqs. (|1^)-(|18|). Due to the spherical symmetry of the band 
we assume Pj{u) = dj{u) = 0. Analogously to the parabolic case, we require that 
eqs. (|17|)"(II1) hold if / depends on k only through e. In the same manner we can 
see that that 

Pij{u) = g{u)5ij , dij{u) = h{u)5ij 



where 



and 



2 uil + aku) 

9W - 



h{u) 



3 m*(l + 2afcM)2 
1 Aak uil + aku) 



m*{l + 2aku) 3m* {l + 2akuY 
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Therefore eq. (0) becomes 

(24) ^ + g{u)V^N - ce|- [g{u)N] + eE h{u)N = -u{u)V . 

Also in this case the equations for and V are equivalent to those obtained by 
means of the spherical harmonics expansion |T^. For = we recover eq. 



6. Homogeneous solutions with constant electric field 

We consider equations (|ll]) and (p^ , and we look for a solution depending only 
on u. Since in this case the electric field must be a constant vector, we choose the 
reference frame so that only the first component of E is different from zero. As a 
consequence the only significant component of V is the first. We limit to assume 
constant kernels % and %o. This fact holds when a homogeneous, intrinsic silicon 
at room temperature is considered. 



It is useful to introduce dimensionless variables. Let be 



w : = 



a 



knTr 



u 
a 



n{w) := u^ilN{u), v{w) := M*£^t*\4(-u) , 



no + 1 



a n '^0 ~ li, 

/3=— -, a:=akU^, ( := cE— 

It is a simple matter that the following equations can be obtained 

dv 

(25) — — = s{w) [an{w + a) + n{w — a)] — [s{w + a) + as {w — a)] n{w) 



(26) C 
where 
(27) 
(28) 

(29) 



-- — {r[w) n) + q{w)n 
aw 



r{w) 
q{w) 



[s{w + a) + as {w — a) + Ps{w)] v{w) . 



- 6{w) ^w{l + aw){l + 2aw) 

2 w{l + aw) 

3 (1 + 2aw;)2 

1 4a w(l + aw) 



1 + 2aw 3 (1 + 2aw)^ 
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As previously mentioned, we recall that we assume the conditions 

n{Q) = v{0) = 0. 

We note that in this case eq. (p6D gives f (0) = as a consequence of the boundary 
condition n{0) = 0. 

We look for solutions of eqs. (p5|)-p6|) satisfying the following conditions 



(30) n(0) = 0, lim v(w) = 

lu— >+oo 

(31) n{w) > for every w > and / n{w) dw > 0. 

Since eqs. (p5D-(p^) are linear and homogeneous, if a solution {n{w),v{w)), satisfying 
the above conditions exists, then, for every c > 0, also {cn{w),cv{w)) is solution. 

It is important to note that the conservation of particle number holds for the 
system (p5D-(p6D. In fact, from eq. ( pSj) it is easy to verify that 

r+oo 

(32) / s{w) [an{w + a) + n{w — a)] — {s{w + a) + as{w — a)) n{w)dw = 0, 
Jo 

for any nonnegative function n such that n{w)s{w) is integrable and n{0) = 0. This 
property forces to choose the boundary on v in (|30D, being t>(0) = 0. 



7. Asymptotic equations 

A simple approximate solution of the previous problem can be obtained in the 
parabolic band case. We introduce a new variable ip defined by 

n{w) = \^ip{w). 

Then, taking into account that now a = 0, eqs. (p5|)-(|26D become 

(33) ~ ay/w + aip{w + a) + \J {w — a)j^il){{w — a)j^ 



— \/w (^w + a + a\J {w — a)+j 



(34) |C Vw^ ^ = ( y^w + a + a^J{w - a)^ + /3y/w )v{w), 

12 



where (z)+ = max{z,0}. 

We look for asymptotic equations of system (p^)-p^) for large values of the energy 
w. To this scope we expand the coefficients of the equations up to the zeroth order 
(for example \Jw + a ~ \/w). We obtain a new set of equations 

(35) - C,'^^ = 'w {[ail)A{w + a) + iI)a{w - a)] - (a + 1)iPa{w)} 

36 va{w) = —— —^w——, 

3(1 + a + p) dw 

where the subscript A indicates the new unknowns. Inserting va{w) in the first 
equation and neglecting a coefficient proportional to w~^, a single equation in ipA is 
obtained 

2(^2 d'^tpA 

(3'') — I 1 I /ON ^ 2 + «V'a(w + a)+ iIja{w - a) - (a + I)V^a(w) = 0. 
3(a + 1 + p) aw^ 

This is a linear difference-differential equation. It is easy to see that exp(Aw) is a 
solution of eq. ( P7| ) if and only if A satisfies the transcendent equation 

(38) _^A^ + ,,-+e---(«+l)=0. 

In Appendix A, we prove that this equation admits only two solutions: A = and 

A G ( — 1,0). We give also a simple approximation for A. 

By means of this solution of eq. (p7|), we get the original quantities 



(39) riAiw) = cV^e'- , ^^(t^) = Ace^" 

3(1 + a + (j) 



where c is a positive constant, in order to satisfy the condition (0). Despite this 
solution is obtained for large values of w it often agrees with the numerical results 
(which we present in next section) in the whole range of w. 



In the limit case C = (no electric field) A = — 1 and the solution of ( p7D gives 
the correct function N obtained from the equilibrium maxwellian distribution func- 
tion lo . 

13 



For high electric fields approximate solutions are known (see 0, for example). 
These are obtained by means of an asymptotic expansion of the original differential- 
difference equations, which are transformed in simple ordinary differential equations 
and are explicitly solved. Our solution, obtained by the asymptotic equation, is 
different with respect to the previous ones, because we do not approximate the 
difference terms /o(e ± hu) {u is the constant optical phonon frequency), appearing 
in the equations, using Taylor formula. In fact the asymptotic equation remains of 
differential- difference type. 

We are interested in calculating the hydrodynamical velocity of electrons 

r+oo 

/ Vi (m) du 

J Un 



N{u) du 





r+oo 

' v{w)dw 






m* J 


/■ + 00 

' n{w)dw 





By a simple calculation, it is possible to verify that, using (|39|), the following value 
is obtained 

_ [k^ 4 f— 

l''^''' - V m* 3(l + a + /5)v/^^Vl^l- 

If we use the approximate value (see Appendix A) of A, we obtain the saturation 

velocity 



r I I - /"j«±ii 8a(a - 1) 
^ iirn^ \VA, I - y ^ ^ 37r(l + a + /?) 



8. Finite energy model 

The difficulties to solve numerically eqs. (p5|)-(|26|) are given by the deviating argu- 
ments in the right side of equations and by the infinite interval for the variable w. 

The first suggests to use a difference scheme to discretize the equations with the 

a 

constraint that the step Aw be such that is an integer. The second difficulty 

Aw 

is solved by considering a finite interval for the variable w. This truncation must 
be made carefully. To explain this aspect of the problem, we recall that we are 
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looking for a solution which describes the no-runaway phenomena. It is a particular 
configuration of the gas, where the applied electric field is balanced by the effect 
of the collisions with the background. In other words, the energy that an electron 
receives from the external electric field is transferred to the phonons so that a spatial 
homogeneous and static configuration is maintained. 

In order to find this solution also in the case of a finite range for the energy w, 
we introduce a cut in the kernels % and %o. Let us choose a positive integer N. Put 
M = Nhuj, we consider the interval [0, M] and the function H{u) which coincides 
with the Heaviside step function 6{u) except for u = where H{u) = 0. Then the 
new kernels are defined by 



These kernels imply that a particle has zero probability of collision in the following 
two cases: 

• the particle before the collision has an energy less than M, but after it should 
have an energy equal or greater than M; 

• the particle before the collision has an energy equal or greater than M. 
Using the new kernels eqs. (p5|)-(|26|) are substituted by the following 



3C^(5(k),e(k)) 



X{e{k),e(k))H{M - max{e(k), e(k)}) , 



Xo{e{k),e(k))H{M - max{£(k), e(k)}) . 



(40) 



-( =A{n){w) 



(41) 



C —-: — {r{w) n) + q{w)n 



dw 



where 



A{n){w) := s{w) aH{aN — w — a)n{w + a) + H{aN — w)n{w — a) 



H{aN — w — a)s{w + a) + aH{aN — w)s{w — a) n{w) 
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B{v){w) := ^H{aN — w — a)s{w + a) + H{aN — w) [as{w — a) + /5s(w)] j v{w) . 
The corresponding boundary conditions becomes 

(42) n(0) = 0, v{aN)=0 



raN 

(43) n{w) > for every w > and / n{w) dw > 0. 

Jo 

The boundary condition on v is the natural consequence of the condition (|30D. The 
choice of the above boundary conditions and the cut are suggested by a simple 
physical picture. No-runaway phenomena may be considered as the limit for t — > 
+00 of a suitable spatial homogeneous solution of eqs. (pT]) , (^Tf). For time- depending 
solutions, the new kernels and the condition v{aN) = imply that a particle with 
initial energy less than M will have, for all time, energy less than M. For the 
same reason particles having energy greater than M will have, for all time, energy 
greater than M. Then, the total number of the particles having energy less than 
M will be constant in time. Different choices can violate the conservation of the 
particle number for electrons having energy less than M. In these cases the number 
of particles having energy less than M could tend to zero or to infinity. 

9. Numerical solutions 

In order to discretize eqs. (140|)- (^If) , we fix the step Aw and consider in the interval 
[0, aN] the points 

Wi = i X Aw (i = 0, 1, 2, 3M") 

aN 

where !N = — — . If we integrate eqs. (pOD-(p]) in the interval [wi, Wi+i] exactly where 



Aw 

it is possible and using trapezoidal rule otherwise, we obtain the linear algebraic 
system 

(44) ^(y.^^-y.)+Ai+,+Ai=0 (z = 0,l,...,X-l), 

Aw; 

(45) vo = 

16 



2C 

(46) -— (rii+iri+i - n^rj) - ({qi+iUi+i + qiUi) - Bi+i - Bi = 

Aw 

(z = l,2,...,K-l), 



where, for any function ip{w), ipi = (p{wi). The corresponding boundary conditions 
are 



(47) no = 0, vj, = 

Aw ^"^ 

(48) > i = 0, and — J2 (ni+i+rii) > 0. 

2 i=o 



Equations (H^) are obtained considering only the interval [wi, Wi+i] for i > 1, because 



in the first point Wq we have written the exact equation t>o = (^e. eq. (^)), which 
derives from eq. ( PD with the boundary condition ?t,(0) = 0. It is simple to deduce 
from (^if ) the following equations 



Aw 

(49) Vj=vy, + — J2{A + A,+,) (j = 0,l,...,X-l). 

Moreover, as a consequence of the appropriate choice of the cut, the conservation of 
the particle number holds (see ref. [§]), being 

(50) 5](A, + A,+i) =0. 

i=0 



Eq. (pOD is the corresponding of eq. (^) introducing the cut in the kernels and 
performing the integral using the trapezoidal rule. 

Now, eq. (^9]) gives Vq = for j = 0. Due to the boundary condition Vj^ = 0, there 
follows that ( ^51) is a consequence of eqs. (^), or equivalently eq. (|49| ) for j = 
coincides with eq. (|45|). For numerical convenience, we prefer to eliminate the case 



j = in (^91) . Now, we must add the condition (H). To select an unique solution, 
we choose the equation ni = 1. This, with the conditions > for i = 2,3, ...!N 
guarantees (^8]). This makes the number of equations equal to the number of the 
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unknowns. For the sake of clarity, we rewrite the complete linear system 



(51) 



V7^ = 0, Vo = 0, no = 0, ni = l, 



(52) = ^^ + ^^(A, + A,+i) (j = l,2,...,N-l) 

i=j 



(53) 



2C 



[ni+irj+i - niVi] - ((g^+i^i+i + qirii) - Bi+i - Bi = 



Aw 

(2 = 1,2,. ..,N-1). 

10. Numerical results and conclusion 

We are interested into solve eqs. (^OD -(^ID in the case of a silicon bulk device. The 
appropriate values for the parameters are given in the following table. 



m* = 


0.32 me 


Tl = 


300 K 


huj = 0.063 eV 


X = 




DtK 


= 11.4 eV 


p = 2330 Kg m-3 


Xq = 




= 


9 eV 


Vq = 9040 m sec~^. 


Oik = 


0.5 eV-^ 









Here, mg denotes the electron rest mass. Using these parameters, we get a ~ 2.437 
and (3 ^ 5.986. The values for the electric field are |E| = 10'^, 10^ and 10^ Vcm^^. 
Except for low electric field there is a great difference between P (parabolic) and NP 
(non-parabolic) band. In fig. 1 we plot the velocity versus the electric field, obtained 
using the asymptotic solutions (^) and values yielded by numerical integration of 
eqs. (^0|)-(41) both in P and NP cases. In the first case the value for the saturation 
velocity obtained is 1.27 x 10^ m/ sec. The second model does not show a limit value 
as in the Monte Carlo simulations [Q. For the largest value of |E| the velocity is 
0.95 X 10^ mj sec The experimental value is 1 x 10^ m/sec. In order to compare 
numerical and asymptotic solutions, we have chosen the condition 



(54) 



aN i-aN 

nA{w) dw = / n{w) dw 
Jo 



where n^(w) indicates the asymptotic solution. The integral in ( |5^ of the function 
Ha is performed analytically and the second numerically. 

In such way the constant c in (^Up is chosen. The meaning of eq. ( |5^ ) is evident and 
coherent with the criterion of the choice of the cut. We have made many numerical 
experiments, varying the values of the electric field, the step and the cut N . In 
fig. 2 the quantity N{u), both for P and NP approximation, is shown for a low value 
of the electric field. Nevertheless a small range of the energy — Ahu is sufficient 
to contain almost all the electrons, the numerical results for P case agree with the 
asymptotic solution. This is due to the fact that for a small values of the electric 
field A ^ — 1. The quantity Vi{u) is shown in fig. 3. Here the asymptotic solution 
is given by the formula 

VAyw) = -—^ r- 

3 Ww + a + a\/w — a + P^/w] 

This choice is not coherent with the expansion, but it allows a better agreement 
with the numerical solutions of P case and it puts in evidence a discontinuity in 
the derivative at the point w = a <^=^ u = hu, due to the term y^w — a. This 
square root is also present in the complete set of equations (|25|)-(p6D. Then, such 
discontinuity is also expected in the solutions. A deeper analysis could be necessary 
to prove mathematically this fact. Figs. 4 and 5 shown an intermediate case where 
the agreement is poor. 

When the electric field is high (figs. 6 and 7) very small variations are observable. 
This agreement was reasonable for large values of u. In the other cases this is due 
to the number of particles having low energy which is enough small with respect to 
that of the whole gas. 

The last figures show a quite surprising result in the parabolic case. We have used 
the same value of the electric field E as in fig. 6 and 7. Here, a small value of N is 
used. The density agrees with the asymptotic solutions. The quantity Vi subject 
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to the boundary condition Vi{Na) — 0, follows the asymptotic solutions almost to 
the end of the interval. 
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Appendix A 

The function 

for C > is convex, because ry"(A) > for every real A. 

Since r){0) — and r]'{0) — a{a — 1) > 0, then there follows that the minimum is 
negative and is reached at a point < 0. 

Prom r]{-l) = 2^ (3(a + 1 + f3)y'^ (recall that a = e") we get -1 < A^ < 0. Then, 
there exists only X ^ such that 77(A) = 0. Moreover A e [— 1, 0[. An approximation 
of this value is obtained expanding e^^" around A = up the term of second order. 
A second order algebraic equations follows 



+ — a+1 



3{a+l + (3) 2 



A^ + a(a- 1)A = 0. 



Then 

6{a + 1 + P)a{a - 1) 



A~ 



4:C + 3{a + l + P)a^{a + l) 
Por C ^ 1 we get 

- 3a(a-l)(a + l + (3) 
A ~ . 

2e 
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Fig. 1 Hydrodynamical velocity (m/sec) versus electric field (V/cm). Solid line: results using asymptotic 
solution. Circles: results using numerical solution in the parabolic case. Disks: results using numerical 
solution in the non parabolic case. 
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Fig. 2 The function N for E = lO^V/ cm. Solid line: asymptotic solution. Dashed line: numerical solution 
- parabolic case -. Circles: numerical solution - Kane case -. 
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Fig. 3 The function Vi for E = lO^V/ cm. Solid line: asymptotic solution. Dashed line: numerical solution 
- parabolic case -. Circles: numerical solution - Kane case -. 
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Fig. 5 The function Vi for E = W^V/cm. Solid line: asymptotic solution. Dashed line: numerical solution 
- parabolic case -. Circles: numerical solution - Kane case -. 
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Fig. 7 The function Vi for E = lO^V/ cm. Solid line: asymptotic solution. Dashed line: numerical solution 
- parabolic case -. Circles: numerical solution - Kane case -. 



25 



